A. Objetivos:
import numpy as np
import cv2 as cv
import pywt
import pywt.data
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
B. Monte o seu google drive e verifique os seus arquivos:
from google.colab import drive
drive.mount('/content/gdrive')
C. Codificação de Luminância (P&B) com DWT para a Pimentas
#img = mpimg.imread('/content/gdrive/My Drive/UFABC2019/CSM/lab4/lucas.jpeg')
img = mpimg.imread('/content/gdrive/My Drive/UFABC2019/CSM/lab4/bruno.png')
# Conversão de BGRA (formato original da imagem) para RGB
img = cv.cvtColor(img, cv.COLOR_BGRA2BGR)
img_gray = cv.cvtColor(img, cv.COLOR_RGB2GRAY)
coefs2 = pywt.dwt2(img_gray,'haar', mode='periodization') #1 nÃvel de DWT
(cA, (cH, cV, cD)) = coefs2 #Separando os coeficientes
imgr = pywt.idwt2(coefs2, 'haar', mode = 'periodization') #1 nÃvel de IDWT
plt.figure(figsize=(10,10))
plt.subplot(2,2,1)
plt.imshow(cA,'gray'); plt.title("CA - Aproximação")
plt.subplot(2,2,2)
plt.imshow(cV,'gray'); plt.title("CV - Bordas Verticais")
plt.subplot(2,2,3)
plt.imshow(cH,'gray'); plt.title("CH - Bordas Horizontais")
plt.subplot(2,2,4)
plt.imshow(cD,'gray'); plt.title("CD - Bordas Diagonais")
C.1 Cálculo do Erro Quadrático Médio (MSE) e da Relação Sinal RuÃdo de Pico (PSNR)
a) A MSE é obtida calculando somando-se o erro quadrático de reconstrução pixel a pixel entre a Imagem Original (O) da ReconstruÃda (R) e normalizando pela dimensão (LxA) da imagem:
$ MSE = \frac{1}{LA}{\sum_{i=0}^L}{\sum_{j=0}^A [O(i,j) - R(i,j)]^2}$
b) A SNR de pico (PSNR) é definida para cada plano componente da imagem como:
$ PSNR = 10.log_{10} \left( \frac{{MAX_I}^2}{MSE} \right) $
sendo $MAX_I$ o valor máximo do pixel, que para 8 bits equivale a 255, logo:
$ PSNR = 20.log_{10}(255) - 10.log_{10} (MSE) $
OBS.: Para uma imagem RGB, $ MSE = MSE_R + MSE_G + MSE_B $, sendo similar definiação para YCrCb e HSV
# Calculo da MSE P&B
A, L, Camadas = img.shape
dif = img_gray - imgr
MSE_gray = np.sum(np.matmul(dif,np.transpose(dif)))/(A*L)
print("MSE_Y = {:.2e}".format(MSE_gray))
PSNR_Y = 20*np.log10(255) - 10*np.log10(MSE_gray)
print("PSNR_Luma = {:.2f} dB".format(PSNR_Y))
plt.figure(figsize=(20,10))
infograf = "Imagem ReconstruÃda de Luminância (Y) com PSNR = " + str(np.uint8(PSNR_Y)) + ' dB'
plt.subplot(1,2,1); plt.imshow(img_gray,'gray'); plt.title("Imagem Original P&B")
plt.subplot(1,2,2); plt.imshow(imgr,'gray'); plt.title(infograf)
D. Teste das Funções de Multiresolução wavedec2() e waverec2()
C = pywt.wavedec2(img_gray,'haar', mode = 'symmetric', level=2) # Dois nÃveis de decomposição DWT
imgr2 = pywt.waverec2(C, 'haar', mode = 'symmetric') # Dois nÃveis de IDWT
# Para extrair os coeficientes de cada nÃvel
cA2 = C[0] # Coeficientes de Aprosimação nÃvel 2
(cH1, cV1, cD1) = C[-1] # Coeficientes de Detalhes nÃvel 1
(cH2, cV2, cD2) = C[-2] # Coeficientes de Detalhes nÃvel 2
# Imagem Original
img_gray = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
# Plot dos coeficientes do nÃvel 2
plt.figure(figsize=(5,5))
plt.subplot(2,2,1)
plt.imshow(cA2, 'gray'); plt.title('Ap. N2: CA2')
plt.subplot(2,2,2)
plt.imshow(cV2, 'gray'); plt.title('B. V. N2: CV2')
plt.subplot(2,2,3)
plt.imshow(cH2, 'gray'); plt.title('B. H. N2: CH2')
plt.subplot(2,2,4)
plt.imshow(cD2, 'gray'); plt.title('B. D. N2: CD2')
# Plot Original e Reconstrução
plt.figure(figsize=(20,10))
plt.subplot(1,2,1); plt.imshow(img_gray,'gray'); plt.title('Imagem Original')
plt.subplot(1,2,2); plt.imshow(imgr2,'gray'); plt.title('Imagem ReconstruÃda')
E. Efetuar uma "Montagem" com wavedec2() e wavedecn()
1º NÃvel
CV1 = cV1.copy()
CH1 = cH1.copy()
CD1 = cD1.copy()
2º NÃvel
CA2 = cA2.copy()
CH2 = cH2.copy()
CV2 = cV2.copy()
CD2 = cD2.copy()
# Matriz Final Completa
CC1 = np.bmat([[CA2,CV2],[CH2,CD2]])
h,w = CV1.shape
CC1 = CC1[:h,:w]
CC = np.bmat([[CC1,CV1],[CH1,CD1]])
print(CA1.shape, CH1.shape, CV1.shape, CD1.shape)
plt.figure(figsize=(20,20))
plt.imshow(CC,'gray')
plt.title('Codificação de Imagem em multinÃvel com função wavedec2()')
F. Reconstrução de Imagem Colorida
# Imagem Original
plt.figure(figsize=(20,20))
plt.imshow(img); plt.title("Imagem Original")
# Codificação por planos de cores
# Plano Vermelho
coefs_R = pywt.dwt2(img[:,:,0],'haar', mode='periodization') #1 nÃvel de DWT R
(cA_R, (cH_R, cV_R, cD_R)) = coefs_R #Separando os coeficientes
cr_R = pywt.idwt2(coefs_R, 'haar', mode = 'periodization') #1 nÃvel de IDWT R
plt.figure(figsize=(20,20))
plt.subplot(2,2,1)
plt.imshow(cA_R,'Reds_r'); plt.title("CA_Red - Aproximação")
plt.subplot(2,2,2)
plt.imshow(cV_R,'Reds_r'); plt.title("CV_Red - Bordas Verticais")
plt.subplot(2,2,3)
plt.imshow(cH_R,'Reds_r'); plt.title("CH_Red - Bordas Horizontais")
plt.subplot(2,2,4)
plt.imshow(cD_R,'Reds_r'); plt.title("CD_Red - Bordas Diagonais")
plt.figure(figsize=(20,20))
plt.imshow(cr_R, 'Reds_r'); plt.title("Imagem ReconstruÃda Red")
# Plano Verde
coefs_G = pywt.dwt2(img[:,:,1],'haar', mode='periodization') #1 nÃvel de DWT G
(cA_G, (cH_G, cV_G, cD_G)) = coefs_G #Separando os coeficientes
cr_G = pywt.idwt2(coefs_G, 'haar', mode = 'periodization') #1 nÃvel de IDWT G
plt.figure(figsize=(20,20))
plt.subplot(2,2,1)
plt.imshow(cA_G,'Greens_r'); plt.title("CA_Green - Aproximação")
plt.subplot(2,2,2)
plt.imshow(cV_G,'Greens_r'); plt.title("CV_Green - Bordas Verticais")
plt.subplot(2,2,3)
plt.imshow(cH_G,'Greens_r'); plt.title("CH_Green - Bordas Horizontais")
plt.subplot(2,2,4)
plt.imshow(cD_G,'Greens_r'); plt.title("CD_Green - Bordas Diagonais")
plt.figure(figsize=(20,20))
plt.imshow(cr_G, 'Greens_r'); plt.title("Imagem ReconstruÃda Green")
# Plano Azul
coefs_B = pywt.dwt2(img[:,:,2],'haar', mode='periodization') #1 nÃvel de DWT B
(cA_B, (cH_B, cV_B, cD_B)) = coefs_B #Separando os coeficientes
cr_B = pywt.idwt2(coefs_B, 'haar', mode = 'periodization') #1 nÃvel de IDWT B
plt.figure(figsize=(20,20))
plt.subplot(2,2,1)
plt.imshow(cA_B,'Blues_r'); plt.title("CA_Blue - Aproximação")
plt.subplot(2,2,2)
plt.imshow(cV_B,'Blues_r'); plt.title("CV_Blue - Bordas Verticais")
plt.subplot(2,2,3)
plt.imshow(cH_B,'Blues_r'); plt.title("CH_Blue - Bordas Horizontais")
plt.subplot(2,2,4)
plt.imshow(cD_B,'Blues_r'); plt.title("CD_Blue - Bordas Diagonais")
plt.figure(figsize=(20,20))
plt.imshow(cr_B, 'Blues_r'); plt.title("Imagem ReconstruÃda Blue")
# Reconstrução NÃvel 1 Colorida
A1, L1 = cA_R.shape
imgrec1 = np.zeros((A1,L1,3))
imgrec1[:,:,0] = cA_R.copy()
imgrec1[:,:,1] = cA_G.copy()
imgrec1[:,:,2] = cA_B.copy()
plt.figure(figsize=(10,10))
plt.imshow(imgrec1); plt.title("Imagem ReconstruÃda DWT/IDWT NÃvel 1")
G. Salvando as Aproximações e depois fazendo download dos arquivos, calcular a taxa de compressão com o original
cv.imwrite('/content/gdrive/My Drive/UFABC2019/CSM/lab4/bruno_DWT_N1_Y.bmp', cA) # Aproximação NÃvel 1 só Y
cv.imwrite('/content/gdrive/My Drive/UFABC2019/CSM/lab4/bruno_DWT_N2_Y.bmp', cA2) # Aproximação NÃvel 2 só Y
H. Gravando o Arquivo Codificado DWT/IDWT nÃvel 1 Colorido, calcular a taxa de compressão com o original
cv.imwrite('/content/gdrive/My Drive/UFABC2019/CSM/lab4/bruno_DWT_N1_colorida.bmp', imgrec1) # Gravando Aproximação NÃvel 1 Colorida
I. Reconstrução da Imagem colorida e Cálculo da MSE de cada plano de cor e da PSNR total
# Reconstrução Colorida Original
A,L = imgr2.shape
imgrec = np.zeros((A,L,3))
imgrec[:,:,0] = cr_R.copy()
imgrec[:,:,1] = cr_G.copy()
imgrec[:,:,2] = cr_B.copy()
# Calculo do MSE colorida
dif2 = img - imgrec
MSE_R = np.sum(np.matmul(dif2[:,:,0],np.transpose(dif2[:,:,0])))/(A*L) # Erro Quadrático Médio plano R
MSE_G = np.sum(np.matmul(dif2[:,:,1],np.transpose(dif2[:,:,1])))/(A*L) # Erro Quadrático Médio plano G
MSE_B = np.sum(np.matmul(dif2[:,:,2],np.transpose(dif2[:,:,2])))/(A*L) # Erro Quadrático Médio plano B
print("MSE_Red= {:.2e}".format(MSE_R), " MSE_Green= {:.2e}".format(MSE_G), " MSE_Blue= {:.2e}".format(MSE_B))
# Cálculo da SNR de pico colorida (PSNR), 3 camadas R, G e B
PSNR = 20*np.log10(255) - 10*np.log10(MSE_R + MSE_G + MSE_B)
print("PSNR total = {:.2f} dB".format(PSNR))
plt.figure(figsize=(20,20))
infograf2 = "Imagem ReconstruÃda com PSNR total = " + str(np.uint8(PSNR)) + ' dB'
plt.imshow(imgrec); plt.title(infograf2)
# Imagem do grupo
img = cv.imread('/content/gdrive/My Drive/UFABC2019/CSM/lab4/grupo.jpg')
# Convertendo para YCrCb e extraindo canais
img_YCrCb = cv.cvtColor(img, cv.COLOR_BGR2YCrCb)
#Y,Cr,Cb = cv.split(img_YCrCb)
Y,Cr,Cb = img_YCrCb[:,:,0], img_YCrCb[:,:,1], img_YCrCb[:,:,2]
print("Componente Y")
Y_DWT = pywt.dwt2(Y,'haar', mode='periodization') #1 nÃvel de DWT
(cA, (cH, cV, cD)) = Y_DWT
Y_IDWT = pywt.idwt2(Y_DWT, 'haar', mode = 'periodization') #1 nÃvel de IDWT
plt.figure(figsize=(18,8))
plt.subplot(2,2,1);plt.imshow(cA,'gray'); plt.title("CA - Aproximação")
plt.subplot(2,2,2);plt.imshow(cV,'gray'); plt.title("CV - Bordas Verticais")
plt.subplot(2,2,3);plt.imshow(cH,'gray'); plt.title("CH - Bordas Horizontais")
plt.subplot(2,2,4);plt.imshow(cD,'gray'); plt.title("CD - Bordas Diagonais")
print("Componente Cb")
Cb_DWT = pywt.dwt2(Cb,'haar', mode='periodization') #1 nÃvel de DWT
(cA, (cH, cV, cD)) = Cb_DWT
Cb_IDWT = pywt.idwt2(Cb_DWT, 'haar', mode = 'periodization') #1 nÃvel de IDWT
plt.figure(figsize=(18,8))
plt.subplot(2,2,1);plt.imshow(cA,'gray'); plt.title("CA - Aproximação")
plt.subplot(2,2,2);plt.imshow(cV,'gray'); plt.title("CV - Bordas Verticais")
plt.subplot(2,2,3);plt.imshow(cH,'gray'); plt.title("CH - Bordas Horizontais")
plt.subplot(2,2,4);plt.imshow(cD,'gray'); plt.title("CD - Bordas Diagonais")
print("Componente Cr")
Cr_DWT = pywt.dwt2(Cr,'haar', mode='periodization') #1 nÃvel de DWT
(cA, (cH, cV, cD)) = Cr_DWT
Cr_IDWT = pywt.idwt2(Cr_DWT, 'haar', mode = 'periodization') #1 nÃvel de IDWT
plt.figure(figsize=(18,8))
plt.subplot(2,2,1);plt.imshow(cA,'gray'); plt.title("CA - Aproximação")
plt.subplot(2,2,2);plt.imshow(cV,'gray'); plt.title("CV - Bordas Verticais")
plt.subplot(2,2,3);plt.imshow(cH,'gray'); plt.title("CH - Bordas Horizontais")
plt.subplot(2,2,4);plt.imshow(cD,'gray'); plt.title("CD - Bordas Diagonais")
# Calculo da MSE P&B
A, L, Camadas = img_YCrCb.shape
imgr = np.zeros(img_YCrCb.shape, dtype=Y_IDWT.dtype)
imgr[:A,:L,0] = Y_IDWT[:A,:L]
imgr[:A,:L,1] = Cr_IDWT[:A,:L]
imgr[:A,:L,2] = Cb_IDWT[:A,:L]
def MSE(A,B):
difference_array = np.subtract(A, B)
squared_array = np.square(difference_array)
mse = squared_array.mean()
return mse
def PSNR(mse):
return 20*np.log10(255) - 10*np.log10(mse)
MSE_YCrCb = MSE(img_YCrCb.astype(np.float64), imgr)
print("MSE_YCrCb = {:.2e}".format(MSE_YCrCb))
PSNR_YCrCb = PSNR(MSE_YCrCb)
print("PSNR_YCrCb = {:.2f} dB".format(PSNR_YCrCb))
plt.figure(figsize=(18,8))
plt.subplot(1,2,1); plt.imshow(img_YCrCb,'gray'); plt.title("Imagem YCrCb")
infograf = "Imagem ReconstruÃda (YCrCb) com PSNR = {:.2f} dB".format(PSNR_YCrCb)
plt.subplot(1,2,2); plt.imshow(np.uint8(imgr),'gray'); plt.title(infograf)